—– ANALYSES ——
Analysis 1: the association between sv2a covariance and functional
connectivity matrix
# generate the sv2s covariance matrix in abeta positive group
w.cov.ad <- cor(as.matrix(df.ad[,paste0("W.SDM8.V",1:200)]),method = "spearman")
w.cov.z.ad <- fisherz(w.cov.ad)
w.cov.z.ad[which(w.cov.z.ad==Inf)] <- NA
corrplot(w.cov.z.ad, diag = FALSE, tl.pos = "n", tl.cex = 0.5, method = "color",is.corr = F)
# generate the sv2s covariance matrix in cognitively normal abeta negative group
w.cov.cn0 <- cor(as.matrix(df.cn0[,paste0("W.SDM8.V",1:200)]),method = "spearman")
w.cov.z.cn0 <- fisherz(w.cov.cn0)
w.cov.z.cn0[which(w.cov.z.cn0==Inf)] <- NA
corrplot(w.cov.z.cn0, diag = FALSE, tl.pos = "n", tl.cex = 0.5, method = "color",is.corr = F)
# prepare data fram for analyses
w.cov.vecterized.ad <- w.cov.z.ad[lower.tri(w.cov.z.ad)]
w.cov.vecterized.cn0 <- w.cov.z.cn0[lower.tri(w.cov.z.cn0)]
fc.verterized <- fc.mat[lower.tri(fc.mat)]
df.stat.cov <- data.frame(w.cov.ad= w.cov.vecterized.ad,
w.cov.cn0=w.cov.vecterized.cn0,
fc=fc.verterized)
## test for Abeta+ group
mod.lm.ad <- lm(data = df.stat.cov,w.cov.ad~fc)
mod.lm.ad <- summary(lm.beta(mod.lm.ad))
beta.val.ad <- mod.lm.ad$coefficients["fc","Standardized"]
mod.lm.ad
##
## Call:
## lm(formula = w.cov.ad ~ fc, data = df.stat.cov)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.46136 -0.06379 0.00019 0.06460 0.41015
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -3.695e-04 NA 6.820e-04 -0.542 0.588
## fc 7.572e-06 7.195e-05 7.461e-04 0.010 0.992
##
## Residual standard error: 0.0962 on 19898 degrees of freedom
## Multiple R-squared: 5.176e-09, Adjusted R-squared: -5.025e-05
## F-statistic: 0.000103 on 1 and 19898 DF, p-value: 0.9919
beta.shuffled.ad <- mapply(function(x){
fc.mat.shuffled <- shuffled.fc[[x]]
df.cov.shuffled <- data.frame(w.cov.ad= w.cov.vecterized.ad,
w.cov.cn0=w.cov.vecterized.cn0,
fc=fc.mat.shuffled[lower.tri(fc.mat.shuffled)])
mod.lm.shuffled <- lm(data = df.cov.shuffled,w.cov.ad~fc)
mod.lm.shuffled <- summary(lm.beta(mod.lm.shuffled))
beta.val.shuffled <- mod.lm.shuffled$coefficients[2,2]
}, 1:1000)
p.val.shuffled.ad <- length(which(beta.shuffled.ad > beta.val.ad))/length(beta.shuffled.ad)
## test for CN abeta- group
mod.lm.cn0 <- lm(data = df.stat.cov,w.cov.cn0~fc)
mod.lm.cn0 <- summary(lm.beta(mod.lm.cn0))
beta.val.cn0 <- mod.lm.cn0$coefficients["fc","Standardized"]
mod.lm.cn0
##
## Call:
## lm(formula = w.cov.cn0 ~ fc, data = df.stat.cov)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.71767 -0.11431 -0.00083 0.11342 0.72274
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) 0.001930 NA 0.001221 1.580 0.114
## fc -0.001896 -0.010061 0.001336 -1.419 0.156
##
## Residual standard error: 0.1723 on 19898 degrees of freedom
## Multiple R-squared: 0.0001012, Adjusted R-squared: 5.098e-05
## F-statistic: 2.014 on 1 and 19898 DF, p-value: 0.1558
beta.shuffled.cn0 <- mapply(function(x){
fc.mat.shuffled <- shuffled.fc[[x]]
df.cov.shuffled <- data.frame(w.cov.ad= w.cov.vecterized.ad,
w.cov.cn0=w.cov.vecterized.cn0,
fc=fc.mat.shuffled[lower.tri(fc.mat.shuffled)])
mod.lm.shuffled <- lm(data = df.cov.shuffled,w.cov.cn0~fc)
mod.lm.shuffled <- summary(lm.beta(mod.lm.shuffled))
beta.val.shuffled <- mod.lm.shuffled$coefficients[2,2]
}, 1:1000)
p.val.shuffled.cn0 <- length(which(beta.shuffled.cn0 > beta.val.cn0))/length(beta.shuffled.cn0)
## comparison between Abeta+ group and CN abeta- group with bootstrapping
### define function
sv2a.cov.z.boot <- function(mat,ind){
sv2a.mat <- as.matrix(mat[ind,paste0("W.SDM8.V",1:200)])
w.cov.boot <- cor(sv2a.mat,method = "spearman")
w.cov.z.boot <- fisherz(w.cov.boot)
w.cov.z.boot[which(w.cov.z.boot==Inf)] <- NA
w.cov.vecterized.boot<- w.cov.z.boot[lower.tri(w.cov.z.boot)]
fc.verterized <- fc.mat[lower.tri(fc.mat)]
df.cov.boot <- data.frame(w.cov.boot= w.cov.vecterized.boot,
fc=fc.verterized)
mod.lm <- lm(data = df.cov.boot,w.cov.boot~fc)
mod.lm <- summary(lm.beta(mod.lm))
beta.val <- mod.lm$coefficients["fc","Standardized"]
return(beta.val)
}
beta.ad.boot <- boot(data=df.ad,statistic=sv2a.cov.z.boot,R=1000)
beta.cn0.boot <- boot(data=df.cn0,statistic=sv2a.cov.z.boot,R=1000)
mod.t <- t.test(beta.ad.boot$t, beta.cn0.boot$t)
mod.t
##
## Welch Two Sample t-test
##
## data: beta.ad.boot$t and beta.cn0.boot$t
## t = 31.072, df = 1998, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.006763386 0.007674673
## sample estimates:
## mean of x mean of y
## -9.940842e-06 -7.228971e-03
Analysis 2: network spreading Model (Shafiei, et al, Biol Psychiatry
2020)
# define functions
get.network.spreading.beta.spin <- function(mean.sv2a.w,surrogate.maps,fc.thre){
mapply(function(x){
mean.w.ad.surrog.curren <- surrogate.maps[,x]
fc.vals <- fc.mat[as.vector(lower.tri(fc.mat))]
fc.thre.val <- fc.vals[order(fc.vals)][length(fc.vals)*fc.thre+1]
fc.w.sv2a.ad <- mapply(function(k){
fc.tmp <- fc.mat
fc.tmp[which(fc.tmp<fc.thre.val)] <- NA
fc.w.tmp <- mean(mean.w.ad.surrog.curren[-k]*fc.tmp[-k,k],na.rm=T)
fc.w.tmp
}, 1:200)
df.stat.spin <- data.frame(fc.wei.w.ad= fc.w.sv2a.ad,
mean.w = mean.sv2a.w)
mod.lm <- lm(data = df.stat.spin,mean.w~fc.wei.w.ad)
mod.lm <- summary(lm.beta(mod.lm))
beta.val.spin <- mod.lm$coefficients[2,2]
return(beta.val.spin)
}, 1:1000)
}
get.network.spreading.beta.shuffled <- function(mean.sv2a.w, fc.thre){
mapply(function(x){
fc.mat.shuffled <- as.matrix(shuffled.fc[[x]])
fc.shuffled.vals <- fc.mat.shuffled[as.vector(lower.tri(fc.mat.shuffled))]
fc.shuffled.thre.val <- fc.shuffled.vals[order(fc.shuffled.vals)][length(fc.shuffled.vals)*fc.thre+1]
fc.w.sv2a.ad.shuffled <- mapply(function(k){
fc.tmp <- fc.mat.shuffled
fc.tmp[which(fc.tmp<fc.shuffled.thre.val)] <- NA
fc.w.tmp <- mean(mean.sv2a.w[-k]*fc.tmp[-k,k],na.rm=T)
fc.w.tmp
}, 1:200)
df.stat.shuffled <- data.frame(fc.wei.w.ad= fc.w.sv2a.ad.shuffled,
mean.w = mean.sv2a.w)
mod.lm.shuffled <- lm(data = df.stat.shuffled,mean.w~fc.wei.w.ad)
mod.lm.shuffled <- summary(lm.beta(mod.lm.shuffled))
beta.val.shuffled <- mod.lm.shuffled$coefficients[2,2]
return(beta.val.shuffled)
}, 1:1000)
}
# analyses
fc.thre.all <- seq(0,1,0.25)[1:4] # set different functional network sparsity
beta.null.all <- list()
beta.empirical.all <- c()
p.spin.all <- c()
p.shuffled.all <- c()
for (i in 1:4) {
## empirical model
fc.thre <- fc.thre.all[i]
fc.vals <- fc.mat[as.vector(lower.tri(fc.mat))]
fc.thre.val <- fc.vals[order(fc.vals)][length(fc.vals)*fc.thre+1]
fc.w.sv2a.tmp <- mapply(function(x){
fc.tmp <- fc.mat
fc.tmp[which(fc.tmp<fc.thre.val)] <- NA
fc.w.tmp <- mean(mean.w.ad[-x]*fc.tmp[-x,x],na.rm=T)
c(fc.w.tmp)
}, 1:200)
df.stat.tmp <- data.frame(fc.wei.w = fc.w.sv2a.tmp,mean.w = mean.w.ad)
mod.lm <- lm(data = df.stat.tmp,mean.w~fc.wei.w)
mod.lm <- summary(lm.beta(mod.lm))
mod.lm
beta.val.empirical <- mod.lm$coefficients["fc.wei.w","Standardized"]
p.val.empirical <- mod.lm$coefficients["fc.wei.w","Pr(>|t|)"]
## generate null distribution of beta values
null.beta.spin <- get.network.spreading.beta.spin(mean.w.ad,mean.w.ad.surrogates, fc.thre)
null.beta.shuffled <- get.network.spreading.beta.shuffled(mean.w.ad, fc.thre)
df.null.beta <- data.frame(null = c(rep("spin",1000),rep("rewired",1000)),
beta = c(null.beta.spin,null.beta.shuffled))
p.val.spin <- length(which(null.beta.spin > beta.val.empirical))/length(null.beta.spin)
p.val.shuffled <- length(which(null.beta.shuffled > beta.val.empirical))/length(null.beta.shuffled)
beta.null.all[[i]] <- df.null.beta
beta.empirical.all[i] <- beta.val.empirical
p.spin.all[i] <- p.val.spin
p.shuffled.all[i] <- p.val.shuffled
}
beta.empirical.all <- data.frame(beta.empirical.all)
beta.empirical.all$x.point <- c(1,2,3,4)
df.beta.null.all <- data.frame(Sparsity=c(rep("100%",2000),rep("75%",2000),rep("50%",2000), rep("25%",2000)),
null.model=rep(c(rep("spin",1000),rep("shuffled",1000)),4),
beta.val= c(beta.null.all[[1]][,2],beta.null.all[[2]][,2],beta.null.all[[3]][,2],beta.null.all[[4]][,2]))
df.beta.null.all$Sparsity <- factor(df.beta.null.all$Sparsity, levels = c("100%","75%","50%","25%"))
y.p.val <- min(df.beta.null.all$beta.val)
significance <- c(p.spin.all[1],p.shuffled.all[1],p.spin.all[2],p.shuffled.all[2],p.spin.all[3],p.shuffled.all[3],p.spin.all[4],p.shuffled.all[4])
significance <- as.data.frame(significance)
significance$null.model <- c(rep(c("spin","shuffled"),4))
significance$Sparsity <- c(rep(c("100%","75%","50%","25%"),each = 2))
significance$y.postion <- y.p.val
ggplot(data=df.beta.null.all)+
aes(x=Sparsity, y=beta.val)+
geom_point(position = position_jitterdodge(0.65),shape=16, size=0.35, alpha=0.55,aes(fill = factor(null.model),color = factor(null.model)))+
geom_boxplot(notch = F,color="black", alpha=0,aes(fill = factor(null.model),color = factor(null.model)),linewidth=0.25)+
ylab("Regression-derived β-value")+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_rect(color= " black", fill = NA),
axis.text.x = element_text(size=12),
axis.text.y = element_text(size=10),
axis.title.y = element_text(size=14))+
scale_color_manual(values=c("#79AF97FF", "#6A6599FF"))+
geom_point(data=beta.empirical.all,aes(x=x.point,y=beta.empirical.all),colour="#F56F5C",size=2)+
geom_text(data=significance,aes(label = ifelse(significance<0.001, "***",ifelse(significance<0.01,"**",ifelse(significance<0.05,"*",""))),
group = null.model,y = y.postion-0.05),
position = position_dodge(width = .75), vjust = 0.5,size = 16 / .pt)
# subject-level analyses
for (i in 1:4) {
fc.thre <- fc.thre.all[i]
print(paste0("The functional network sparsity at ", fc.thre))
fc.vals <- fc.mat[as.vector(lower.tri(fc.mat))]
fc.thre.val <- fc.vals[order(fc.vals)][length(fc.vals)*fc.thre+1]
beta.network.spreading.sub <- ddply(df.ad,.(SubID),function(x){
w.ad.sub <- x[,paste0("W.SDM8.V",1:200)]
w.ad.sub <- as.vector(as.matrix(w.ad.sub))
fc.w.sv2a.sub <- mapply(function(x){
fc.tmp <- fc.mat
fc.tmp[which(fc.tmp<fc.thre.val)] <- NA
fc.w.tmp <- mean(w.ad.sub[-x]*fc.tmp[-x,x],na.rm=T)
fc.w.tmp
}, 1:200)
df.stat.ad.sub <- data.frame(fc.wei.w.sub = fc.w.sv2a.sub,mean.w = w.ad.sub)
mod.lm <- lm(data = df.stat.ad.sub,mean.w~fc.wei.w.sub)
mod.lm <- summary(lm.beta(mod.lm))
beta.val.sub <- mod.lm$coefficients[2,2]
beta.val.sub
})
beta.network.spreading.sub$variable <- "beta"
t.mod <- t.test(beta.network.spreading.sub$V1)
print(t.mod)
}
## [1] "The functional network sparsity at 0"
##
## One Sample t-test
##
## data: beta.network.spreading.sub$V1
## t = 0.71183, df = 112, p-value = 0.478
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.00832983 0.01767095
## sample estimates:
## mean of x
## 0.004670561
##
## [1] "The functional network sparsity at 0.25"
##
## One Sample t-test
##
## data: beta.network.spreading.sub$V1
## t = -4.3518, df = 112, p-value = 2.996e-05
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.03922860 -0.01468265
## sample estimates:
## mean of x
## -0.02695563
##
## [1] "The functional network sparsity at 0.5"
##
## One Sample t-test
##
## data: beta.network.spreading.sub$V1
## t = -6.8285, df = 112, p-value = 4.66e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.05591255 -0.03076263
## sample estimates:
## mean of x
## -0.04333759
##
## [1] "The functional network sparsity at 0.75"
##
## One Sample t-test
##
## data: beta.network.spreading.sub$V1
## t = -6.3251, df = 112, p-value = 5.346e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.05030069 -0.02630383
## sample estimates:
## mean of x
## -0.03830226
Analysis 3: epicenter connectivity-based sv2a propagation in pooled
abeta+ subjects
# define functions
get.epicenter.by.rank <- function(mean.sv2a.w){
rank.mean.w <- as.vector(rank(mean.sv2a.w))
fc.wei.w <- mapply(function(k){
fc.w.tmp <- mean(mean.sv2a.w[-k]*fc.mat[-k,k])
fc.w.tmp
}, 1:200)
rank.fc.wei.w <- rank(fc.wei.w)
mean.rank.current <- (rank.mean.w+rank.fc.wei.w)/2
## generated null distribution of mean ranks
surrogate.map <- nmnulls$vasa(mean.sv2a.w,atlas='fsLR',density='32k',parcellation=parc,seed=as.integer(45))
mean.rank.null <- mapply(function(j){
mean.w.surrog <- surrogate.map[,j]
rank.mean.w.surrog <- as.vector(rank(mean.w.surrog))
fc.wei.w.surrog<- mapply(function(k){
fc.w.tmp <- mean(mean.w.surrog[-k]*fc.mat[-k,k])
fc.w.tmp
}, 1:200)
rank.fc.wei.w.surrog <- rank(fc.wei.w.surrog)
mean.rank <- (rank.mean.w.surrog+rank.fc.wei.w.surrog)/2
}, 1:1000)
rank.p.vals <- mapply(function(roi){
p.roi <- length(which(mean.rank.current[roi] > mean.rank.null[roi,]))/length(mean.rank.null[roi,])
}, 1:200)
epicenter.roi <- rep(0,n.ROIs)
epicenter.roi[which(rank.p.vals<0.05)] <- 1
epicenter.roi
}
get.epi.pred.beta.shuffled <- function(epicenter.rois, mean.sv2a.w){
mapply(function(j){
fc.mat.shuffled <- shuffled.fc[[j]]
epi.fc.mat.shuffled <- as.matrix(fc.mat.shuffled[epicenter.rois,])
epi.fc.mat.shuffled[which(epi.fc.mat.shuffled==Inf)] <- NA
mean.epi.fc.mat.shuffled <- colMeans(epi.fc.mat.shuffled,na.rm = T)
df.test.ad.shuffled <- data.frame(epi.fc = mean.epi.fc.mat.shuffled,
sv2a.wscore=mean.sv2a.w)
mod.lm.shuffled <- lm(data = df.test.ad.shuffled,sv2a.wscore~epi.fc)
mod.lm.shuffled <- summary(lm.beta(mod.lm.shuffled))
beta.val <- mod.lm.shuffled$coefficients[2,2]
beta.val
}, 1:1000)
}
# Epicenter connectivity-based prediciton of group-level sv2a w score in pooled abeta+ subjects
## get group-level epicenter ROIs
epicenter.ad <- get.epicenter.by.rank(mean.w.ad)
epicenter.ad <- which(epicenter.ad==1)
## epicenter-based prediction
epi.fc.mat.ad <- fc.mat[epicenter.ad,]
mean.epi.fc.ad <- colMeans(epi.fc.mat.ad,na.rm = T)
## prediction
beta.shuffled.ad <- get.epi.pred.beta.shuffled(epicenter.ad,mean.w.ad)
df.stat.epicenter.fc.prediction <- data.frame(epicenter.fc = mean.epi.fc.ad, mean.w = mean.w.ad)
mod.lm <- lm(data = df.stat.epicenter.fc.prediction, mean.w~epicenter.fc)
mod.lm <- summary(lm.beta(mod.lm))
mod.lm
##
## Call:
## lm(formula = mean.w ~ epicenter.fc, data = df.stat.epicenter.fc.prediction)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.35163 -0.09284 0.00231 0.07666 0.31812
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -0.0001211 NA 0.0082469 -0.015 0.988
## epicenter.fc -0.0250385 -0.0617498 0.0287614 -0.871 0.385
##
## Residual standard error: 0.1166 on 198 degrees of freedom
## Multiple R-squared: 0.003813, Adjusted R-squared: -0.001218
## F-statistic: 0.7579 on 1 and 198 DF, p-value: 0.385
beta.val <- mod.lm$coefficients["epicenter.fc","Standardized"]
p.val.shuffled <- length(which(beta.shuffled.ad < beta.val))/length(beta.shuffled.ad)
## plot
x.range <- range(df.stat.epicenter.fc.prediction$epicenter.fc)
y.range <- range(df.stat.epicenter.fc.prediction$mean.w)
tibble(Predicted = df.stat.epicenter.fc.prediction$epicenter.fc, Response = df.stat.epicenter.fc.prediction$mean.w, nets = as.factor(sf.assign$V1)) %>%
ggplot(aes(x=Predicted, y = Response))+
geom_point(shape=21, size=2,aes(fill=factor(nets)))+
geom_smooth(method = "lm",color = "grey55")+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_rect(color= " black", fill = NA),
aspect.ratio = 1)+
xlab("Functional connectivity to epicenter")+
ylab("SV2A-PET W-score")+
scale_fill_manual(values=colors.Schaefer.networks,labels=c("Visual","Motor","DAN","VAN","Limbic","FPCN","DMN"),name="")+
scale_y_continuous(labels = scales::number_format(accuracy = 0.01))+
annotate("text", x = min(x.range)+2/5*diff(x.range), y = max(y.range)+1/8*diff(y.range),
label = paste0(paste0("italic(β)==", round(beta.val,3)), "*\',\'",
"~", expression(italic(p)[italic(rewired)]),ifelse(p.val.shuffled>0.001,paste0("==",round(p.val.shuffled,3)),"<0.001")),
color = "black",parse=T, size=5)
## `geom_smooth()` using formula = 'y ~ x'
# Epicenter connectivity-based prediciton of subject-level sv2a w score in pooled abeta+ subjects
## define subject-specific epicenter
df.epicenters.subject <- ddply(df.ad,.(SubID), function(x){
mean.w.sub <- as.vector(as.matrix(x[,paste0("W.SDM8.V",1:200)]))
epicenters.sub <- get.epicenter.by.rank(mean.w.sub)
epicenters.sub
})
colnames(df.epicenters.subject)[2:201] <- paste0("Epicenter.bin.V",1:200)
df.ad.epicenter <- merge(df.ad,df.epicenters.subject, by ="SubID")
## subject-level prediction
df.stat.sub <- ddply(df.ad.epicenter,.(SubID), function(x){
mean.w.sub <- as.vector(as.matrix(x[,paste0("W.SDM8.V",1:200)]))
epicenters.sub <- which(x[,paste0("Epicenter.bin.V",1:200)]==1)
epi.fc.mat.sub <- fc.mat[epicenters.sub,]
mean.epi.fc.sub <- colMeans(epi.fc.mat.sub,na.rm = T)
fc.quantile.thre <- quantile(mean.epi.fc.sub)
q1.rois <- which(mean.epi.fc.sub<fc.quantile.thre[2])
q2.rois <- which(mean.epi.fc.sub<fc.quantile.thre[3]&mean.epi.fc.sub>=fc.quantile.thre[2])
q3.rois <- which(mean.epi.fc.sub<fc.quantile.thre[4]&mean.epi.fc.sub>=fc.quantile.thre[3])
q4.rois <- which(mean.epi.fc.sub>=fc.quantile.thre[4])
fc.q1 <- mean(mean.epi.fc.sub[q1.rois])
fc.q2 <- mean(mean.epi.fc.sub[q2.rois])
fc.q3 <- mean(mean.epi.fc.sub[q3.rois])
fc.q4 <- mean(mean.epi.fc.sub[q4.rois])
w.q1 <- mean(mean.w.sub[q1.rois])
w.q2 <- mean(mean.w.sub[q2.rois])
w.q3 <- mean(mean.w.sub[q3.rois])
w.q4 <- mean(mean.w.sub[q4.rois])
epi.eu.mat.sub <- eu.d[epicenters.sub,]
mean.epi.eu.d.sub <- colMeans(epi.eu.mat.sub,na.rm = T)
eu.q1 <- mean(mean.epi.eu.d.sub[q1.rois])
eu.q2 <- mean(mean.epi.eu.d.sub[q2.rois])
eu.q3 <- mean(mean.epi.eu.d.sub[q3.rois])
eu.q4 <- mean(mean.epi.eu.d.sub[q4.rois])
df.fc.sub <- data.frame(Quantile = c("Q1","Q2","Q3","Q4"),
Quantile.num = 1:4,
FC.quantile = c(fc.q1,fc.q2,fc.q3,fc.q4),
mean.sv2a = c(w.q1,w.q2,w.q3,w.q4),
EU.dist = c(eu.q1,eu.q2,eu.q3,eu.q4))
})
df.stat.sub <- merge(df.ad.epicenter,df.stat.sub, by ="SubID")
mod.lmer <- lmer(data = df.stat.sub,mean.sv2a ~ FC.quantile + SEX + Age.SDM8 + DX.AB +EU.dist +(1|SubID)) ;summary(mod.lmer)
## boundary (singular) fit: see help('isSingular')
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean.sv2a ~ FC.quantile + SEX + Age.SDM8 + DX.AB + EU.dist +
## (1 | SubID)
## Data: df.stat.sub
##
## REML criterion at convergence: -329
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0007 -0.6760 -0.0116 0.6323 3.0220
##
## Random effects:
## Groups Name Variance Std.Dev.
## SubID (Intercept) 0.00000 0.0000
## Residual 0.02575 0.1605
## Number of obs: 452, groups: SubID, 113
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -3.515e-01 2.815e-01 4.450e+02 -1.249 0.2124
## FC.quantile 2.079e-02 2.761e-02 4.450e+02 0.753 0.4520
## SEXmale 2.864e-02 1.562e-02 4.450e+02 1.833 0.0674 .
## Age.SDM8 -5.009e-04 8.678e-04 4.450e+02 -0.577 0.5641
## DX.ABDementia_1 -1.037e-03 1.799e-02 4.450e+02 -0.058 0.9541
## DX.ABMCI_1 -1.840e-02 1.917e-02 4.450e+02 -0.960 0.3375
## EU.dist 5.086e-02 3.761e-02 4.450e+02 1.352 0.1770
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) FC.qnt SEXmal A.SDM8 DX.ABD DX.ABM
## FC.quantile -0.057
## SEXmale -0.002 -0.010
## Age.SDM8 -0.121 -0.020 0.072
## DX.ABDmnt_1 0.062 0.003 0.041 -0.067
## DX.ABMCI_1 0.040 -0.003 -0.187 0.053 0.428
## EU.dist -0.972 0.062 -0.038 -0.111 -0.077 -0.076
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## plot
y.range <- range(df.stat.sub$mean.sv2a)
ggplot(data= df.stat.sub)+
aes(x=Quantile, y = mean.sv2a)+
geom_point( size=2, color = "#DF8F44FF", alpha=0.5,aes(group=SubID))+
geom_line(color = "#DF8F44FF", alpha=0.5,aes(group=SubID))+
geom_smooth(method = "lm",color = "grey55",aes(x=Quantile.num, y = mean.sv2a))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_rect(color= " black", fill = NA),
axis.text.x = element_text(size=10),
axis.text.y = element_text(size=10))+
xlab(str_wrap("Functional connectivity quantile to subject-specific epicenter", width = 40))+
ylab(str_wrap("SV2A-PET W-score", width = 30))
## `geom_smooth()` using formula = 'y ~ x'
Analysis 4: epicenter connectivity-based sv2a propagation in
subtypes of abeta+ subjects
# hierachical clustering to determine three sv2a w-score subtypes
df.cluster.ad <- df.ad.epicenter
df.cluster.ad.w <- df.cluster.ad[,paste0("W.SDM8.V",1:200)]
df.cluster.ad.w <- scale(t(df.cluster.ad.w))
df.cluster.ad.w <- t(df.cluster.ad.w)
d <- dist(df.cluster.ad.w,method = 'euclidean')
fit.hc <- hclust(d, method = "ward.D")
clusters <- cutree(fit.hc,k=3)
df.ad.epicenter$clusters <- clusters
# group-level epicenter connectivity-based prediction of sv2a w-score for each subtypese
for (i in 1:3) {
print(paste0("Analysis for cluster ", i))
df.ad.epicenter.cluster <- subset(df.ad.epicenter, clusters==i)
mean.w.ad.cluster <- colMeans(df.ad.epicenter.cluster[,paste0("W.SDM8.V",1:200)])
epicenter.ad.cluster <- get.epicenter.by.rank(mean.w.ad.cluster)
epicenter.ad.cluster <- which(epicenter.ad.cluster==1)
## epicenter-based prediction
epi.fc.mat.ad.cluster <- fc.mat[epicenter.ad.cluster,]
mean.epi.fc.ad.cluster <- colMeans(epi.fc.mat.ad.cluster,na.rm = T)
## prediction
beta.shuffled.ad.cluster <- get.epi.pred.beta.shuffled(epicenter.ad.cluster, mean.w.ad.cluster)
df.stat.cluster.epicenter.fc.prediction <- data.frame(epicenter.fc = mean.epi.fc.ad.cluster, mean.w = mean.w.ad.cluster)
mod.lm <- lm(data = df.stat.cluster.epicenter.fc.prediction, mean.w~epicenter.fc)
mod.lm <- summary(lm.beta(mod.lm))
print(mod.lm)
beta.val <- mod.lm$coefficients["epicenter.fc","Standardized"]
p.val.shuffled <- length(which(beta.shuffled.ad.cluster < beta.val))/length(beta.shuffled.ad.cluster)
## plot
x.range <- range(df.stat.cluster.epicenter.fc.prediction$epicenter.fc)
y.range <- range(df.stat.cluster.epicenter.fc.prediction$mean.w)
tibble(Predicted = df.stat.cluster.epicenter.fc.prediction$epicenter.fc, Response = df.stat.cluster.epicenter.fc.prediction$mean.w, nets = as.factor(sf.assign$V1)) %>%
ggplot(aes(x=Predicted, y = Response))+
geom_point(shape=21, size=2,aes(fill=factor(nets)))+
geom_smooth(method = "lm",color = "grey55")+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_rect(color= " black", fill = NA),
aspect.ratio = 1)+
xlab("Functional connectivity to epicenter")+
ylab("SV2A-PET W-score")+
ggtitle(paste0("Epicenter connectivity-based prediction of SV2A for subtype ",i))+
scale_fill_manual(values=colors.Schaefer.networks,labels=c("Visual","Motor","DAN","VAN","Limbic","FPCN","DMN"),name="")+
scale_y_continuous(labels = scales::number_format(accuracy = 0.01))+
annotate("text", x = min(x.range)+2/5*diff(x.range), y = max(y.range)+1/8*diff(y.range),
label = paste0(paste0("italic(β)==", round(beta.val,3)), "*\',\'",
"~", expression(italic(p)[italic(rewired)]),ifelse(p.val.shuffled>0.001,paste0("==",round(p.val.shuffled,3)),"<0.001")),
color = "black",parse=T, size=5)
## sliding window analyses
df.ad.epicenter.cluster$mean.sv2a <- rowMeans(df.ad.epicenter.cluster[,paste0("W.SDM8.V",1:200)])
sub.order.cluster <- order(df.ad.epicenter.cluster$mean.sv2a,decreasing = T)
df.ad.epicenter.cluster.slide <- df.ad.epicenter.cluster[sub.order.cluster,]
num.sub <- length(sub.order.cluster)
win.wid <- 10 # set the width of windows at 10 subjects
win.num <- 9 # set the number of windows at 10 windows
sep <- win.wid-((win.wid*win.num-num.sub)/(win.num-1))
windows.assign <- rollapply(1:num.sub, win.wid, by = sep, c)
for (i.win in 1:win.num) {
df.current.window <- df.ad.epicenter.cluster.slide[as.vector(windows.assign[i.win,]),]
mean.w.current.window <- colMeans(df.current.window[,paste0("W.SDM8.V",1:200)])
epicenters.current.window <- epicenter.ad.cluster
epi.fc.current.window <- fc.mat[epicenters.current.window,]
df.test.current.window <- data.frame(epicenter.fc = colMeans(epi.fc.current.window,na.rm = T),
mean.w = mean.w.current.window)
beta.shuffled.current.window <- get.epi.pred.beta.shuffled(epicenters.current.window,mean.w.current.window)
mod.lm <- lm(data = df.test.current.window,mean.w~epicenter.fc)
mod.lm <- summary(lm.beta(mod.lm))
beta.val <- mod.lm$coefficients["epicenter.fc","Standardized"]
p.val.shuffled <- length(which(beta.shuffled.current.window< beta.val))/length(beta.shuffled.current.window)
r2 <- mod.lm$adj.r.squared
print(paste0("Sliding window analyses for cluster ", i , " window ", i.win))
print(paste0("β=", round(beta.val,3) , ", p.rewired=", p.val.shuffled,", r2=",round(r2,3)))
}
}
## [1] "Analysis for cluster 1"
##
## Call:
## lm(formula = mean.w ~ epicenter.fc, data = df.stat.cluster.epicenter.fc.prediction)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.60106 -0.11353 0.00104 0.14385 0.42430
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) 0.007844 NA 0.014211 0.552 0.582
## epicenter.fc 0.091785 0.115261 0.056215 1.633 0.104
##
## Residual standard error: 0.201 on 198 degrees of freedom
## Multiple R-squared: 0.01329, Adjusted R-squared: 0.008302
## F-statistic: 2.666 on 1 and 198 DF, p-value: 0.1041
##
## [1] "Sliding window analyses for cluster 1 window 1"
## [1] "β=0.051, p.rewired=0.73, r2=-0.002"
## [1] "Sliding window analyses for cluster 1 window 2"
## [1] "β=0.071, p.rewired=0.828, r2=0"
## [1] "Sliding window analyses for cluster 1 window 3"
## [1] "β=0.036, p.rewired=0.698, r2=-0.004"
## [1] "Sliding window analyses for cluster 1 window 4"
## [1] "β=0.073, p.rewired=0.849, r2=0"
## [1] "Sliding window analyses for cluster 1 window 5"
## [1] "β=0.055, p.rewired=0.776, r2=-0.002"
## [1] "Sliding window analyses for cluster 1 window 6"
## [1] "β=0.102, p.rewired=0.933, r2=0.005"
## [1] "Sliding window analyses for cluster 1 window 7"
## [1] "β=0.125, p.rewired=0.959, r2=0.011"
## [1] "Sliding window analyses for cluster 1 window 8"
## [1] "β=0.096, p.rewired=0.912, r2=0.004"
## [1] "Sliding window analyses for cluster 1 window 9"
## [1] "β=0.023, p.rewired=0.619, r2=-0.005"
## [1] "Analysis for cluster 2"
##
## Call:
## lm(formula = mean.w ~ epicenter.fc, data = df.stat.cluster.epicenter.fc.prediction)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.06040 -0.24926 -0.02737 0.23933 0.94678
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) 0.004828 NA 0.025959 0.186 0.853
## epicenter.fc 0.116755 0.077688 0.106482 1.096 0.274
##
## Residual standard error: 0.3662 on 198 degrees of freedom
## Multiple R-squared: 0.006035, Adjusted R-squared: 0.001015
## F-statistic: 1.202 on 1 and 198 DF, p-value: 0.2742
##
## [1] "Sliding window analyses for cluster 2 window 1"
## [1] "β=0.055, p.rewired=0.785, r2=-0.002"
## [1] "Sliding window analyses for cluster 2 window 2"
## [1] "β=0.068, p.rewired=0.831, r2=0"
## [1] "Sliding window analyses for cluster 2 window 3"
## [1] "β=0.047, p.rewired=0.733, r2=-0.003"
## [1] "Sliding window analyses for cluster 2 window 4"
## [1] "β=0.07, p.rewired=0.817, r2=0"
## [1] "Sliding window analyses for cluster 2 window 5"
## [1] "β=0.067, p.rewired=0.827, r2=-0.001"
## [1] "Sliding window analyses for cluster 2 window 6"
## [1] "β=0.024, p.rewired=0.603, r2=-0.004"
## [1] "Sliding window analyses for cluster 2 window 7"
## [1] "β=0.052, p.rewired=0.748, r2=-0.002"
## [1] "Sliding window analyses for cluster 2 window 8"
## [1] "β=0.042, p.rewired=0.687, r2=-0.003"
## [1] "Sliding window analyses for cluster 2 window 9"
## [1] "β=0.07, p.rewired=0.829, r2=0"
## [1] "Analysis for cluster 3"
##
## Call:
## lm(formula = mean.w ~ epicenter.fc, data = df.stat.cluster.epicenter.fc.prediction)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.7511 -0.1424 -0.0174 0.1561 0.5649
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -0.01089 NA 0.01549 -0.703 0.483
## epicenter.fc 0.00197 0.00277 0.05054 0.039 0.969
##
## Residual standard error: 0.2189 on 198 degrees of freedom
## Multiple R-squared: 7.673e-06, Adjusted R-squared: -0.005043
## F-statistic: 0.001519 on 1 and 198 DF, p-value: 0.9689
##
## [1] "Sliding window analyses for cluster 3 window 1"
## [1] "β=-0.054, p.rewired=0.211, r2=-0.002"
## [1] "Sliding window analyses for cluster 3 window 2"
## [1] "β=-0.09, p.rewired=0.103, r2=0.003"
## [1] "Sliding window analyses for cluster 3 window 3"
## [1] "β=-0.076, p.rewired=0.144, r2=0.001"
## [1] "Sliding window analyses for cluster 3 window 4"
## [1] "β=0.063, p.rewired=0.826, r2=-0.001"
## [1] "Sliding window analyses for cluster 3 window 5"
## [1] "β=0.124, p.rewired=0.965, r2=0.01"
## [1] "Sliding window analyses for cluster 3 window 6"
## [1] "β=0.118, p.rewired=0.951, r2=0.009"
## [1] "Sliding window analyses for cluster 3 window 7"
## [1] "β=0.057, p.rewired=0.778, r2=-0.002"
## [1] "Sliding window analyses for cluster 3 window 8"
## [1] "β=-0.071, p.rewired=0.162, r2=0"
## [1] "Sliding window analyses for cluster 3 window 9"
## [1] "β=-0.063, p.rewired=0.175, r2=-0.001"
# subject-level epicenter connectivity-based prediction of sv2a w-score for each subtypese
for (i in 1:3) {
print(paste0("Subject-level analysis for cluster ", i))
df.ad.epicenter.cluster <- subset(df.ad.epicenter, clusters==i)
df.stat.sub.cluster <- ddply(df.ad.epicenter.cluster,.(SubID), function(x){
mean.w.sub <- as.vector(as.matrix(x[,paste0("W.SDM8.V",1:200)]))
epicenters.sub <- which(x[,paste0("Epicenter.bin.V",1:200)]==1)
epi.fc.mat.sub <- fc.mat[epicenters.sub,]
mean.epi.fc.sub <- colMeans(epi.fc.mat.sub,na.rm = T)
fc.quantile.thre <- quantile(mean.epi.fc.sub)
q1.rois <- which(mean.epi.fc.sub<fc.quantile.thre[2])
q2.rois <- which(mean.epi.fc.sub<fc.quantile.thre[3]&mean.epi.fc.sub>=fc.quantile.thre[2])
q3.rois <- which(mean.epi.fc.sub<fc.quantile.thre[4]&mean.epi.fc.sub>=fc.quantile.thre[3])
q4.rois <- which(mean.epi.fc.sub>=fc.quantile.thre[4])
fc.q1 <- mean(mean.epi.fc.sub[q1.rois])
fc.q2 <- mean(mean.epi.fc.sub[q2.rois])
fc.q3 <- mean(mean.epi.fc.sub[q3.rois])
fc.q4 <- mean(mean.epi.fc.sub[q4.rois])
w.q1 <- mean(mean.w.sub[q1.rois])
w.q2 <- mean(mean.w.sub[q2.rois])
w.q3 <- mean(mean.w.sub[q3.rois])
w.q4 <- mean(mean.w.sub[q4.rois])
epi.eu.mat.sub <- eu.d[epicenters.sub,]
mean.epi.eu.d.sub <- colMeans(epi.eu.mat.sub,na.rm = T)
eu.q1 <- mean(mean.epi.eu.d.sub[q1.rois])
eu.q2 <- mean(mean.epi.eu.d.sub[q2.rois])
eu.q3 <- mean(mean.epi.eu.d.sub[q3.rois])
eu.q4 <- mean(mean.epi.eu.d.sub[q4.rois])
df.fc.sub <- data.frame(Quantile = c("Q1","Q2","Q3","Q4"),
Quantile.num = 1:4,
FC.quantile = c(fc.q1,fc.q2,fc.q3,fc.q4),
mean.sv2a = c(w.q1,w.q2,w.q3,w.q4),
EU.dist = c(eu.q1,eu.q2,eu.q3,eu.q4))
})
df.stat.sub.cluster <- merge(df.ad.epicenter.cluster,df.stat.sub.cluster, by ="SubID")
mod.lmer <- lmer(data = df.stat.sub.cluster,mean.sv2a ~ FC.quantile + SEX + Age.SDM8 + DX.AB +EU.dist+(1|SubID)) ;print(summary(mod.lmer))
# plot
y.range <- range(df.stat.sub.cluster$mean.sv2a)
p <- ggplot(data= df.stat.sub.cluster)+
aes(x=Quantile, y = mean.sv2a)+
geom_point( size=2, color = "#DF8F44FF", alpha=0.5,aes(group=SubID))+
geom_line(color = "#DF8F44FF", alpha=0.5,aes(group=SubID))+
geom_smooth(method = "lm",color = "grey55",aes(x=Quantile.num, y = mean.sv2a))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_rect(color= " black", fill = NA),
axis.text.x = element_text(size=10),
axis.text.y = element_text(size=10))+
xlab(str_wrap("Functional connectivity quantile to subject-specific epicenter", width = 40))+
ylab(str_wrap("SV2A-PET W-score", width = 30))+
ggtitle(paste0("Subject-level epicenter connectivity-based prediciton for subtype ",i))
print(p)
}
## [1] "Subject-level analysis for cluster 1"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean.sv2a ~ FC.quantile + SEX + Age.SDM8 + DX.AB + EU.dist +
## (1 | SubID)
## Data: df.stat.sub.cluster
##
## REML criterion at convergence: -134.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.50254 -0.71621 -0.06735 0.63409 2.83313
##
## Random effects:
## Groups Name Variance Std.Dev.
## SubID (Intercept) 0.001155 0.03399
## Residual 0.024739 0.15729
## Number of obs: 208, groups: SubID, 52
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.522823 0.418079 200.667280 -1.251 0.213
## FC.quantile 0.040789 0.040758 155.784056 1.001 0.318
## SEXmale 0.026953 0.025238 47.039425 1.068 0.291
## Age.SDM8 -0.001613 0.001417 47.199081 -1.139 0.261
## DX.ABDementia_1 0.006228 0.027505 47.619268 0.226 0.822
## DX.ABMCI_1 -0.015792 0.035193 47.204933 -0.449 0.656
## EU.dist 0.085504 0.054673 194.805403 1.564 0.119
##
## Correlation of Fixed Effects:
## (Intr) FC.qnt SEXmal A.SDM8 DX.ABD DX.ABM
## FC.quantile -0.052
## SEXmale 0.027 -0.008
## Age.SDM8 -0.207 -0.012 -0.056
## DX.ABDmnt_1 0.057 0.011 -0.241 0.023
## DX.ABMCI_1 -0.042 0.003 -0.258 0.299 0.399
## EU.dist -0.965 0.056 -0.029 -0.053 -0.089 -0.054
## `geom_smooth()` using formula = 'y ~ x'
## [1] "Subject-level analysis for cluster 2"
## boundary (singular) fit: see help('isSingular')
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean.sv2a ~ FC.quantile + SEX + Age.SDM8 + DX.AB + EU.dist +
## (1 | SubID)
## Data: df.stat.sub.cluster
##
## REML criterion at convergence: -25.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7761 -0.5048 -0.1646 0.6790 2.3266
##
## Random effects:
## Groups Name Variance Std.Dev.
## SubID (Intercept) 0.00000 0.000
## Residual 0.02788 0.167
## Number of obs: 72, groups: SubID, 18
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.558297 0.812288 65.000000 -0.687 0.494
## FC.quantile -0.034249 0.070572 65.000000 -0.485 0.629
## SEXmale 0.045913 0.044589 65.000000 1.030 0.307
## Age.SDM8 0.003379 0.003085 65.000000 1.095 0.277
## DX.ABDementia_1 -0.017758 0.050562 65.000000 -0.351 0.727
## DX.ABMCI_1 0.017808 0.052323 65.000000 0.340 0.735
## EU.dist 0.040302 0.104171 65.000000 0.387 0.700
##
## Correlation of Fixed Effects:
## (Intr) FC.qnt SEXmal A.SDM8 DX.ABD DX.ABM
## FC.quantile 0.030
## SEXmale -0.108 -0.027
## Age.SDM8 -0.298 -0.006 0.265
## DX.ABDmnt_1 0.140 -0.015 0.233 -0.221
## DX.ABMCI_1 0.127 0.011 -0.118 0.138 0.234
## EU.dist -0.961 -0.028 0.008 0.026 -0.113 -0.187
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## `geom_smooth()` using formula = 'y ~ x'
## [1] "Subject-level analysis for cluster 3"
## boundary (singular) fit: see help('isSingular')
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean.sv2a ~ FC.quantile + SEX + Age.SDM8 + DX.AB + EU.dist +
## (1 | SubID)
## Data: df.stat.sub.cluster
##
## REML criterion at convergence: -104
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.01054 -0.65480 0.06666 0.55215 2.93787
##
## Random effects:
## Groups Name Variance Std.Dev.
## SubID (Intercept) 0.00000 0.0000
## Residual 0.02605 0.1614
## Number of obs: 172, groups: SubID, 43
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -6.747e-02 4.523e-01 1.650e+02 -0.149 0.882
## FC.quantile 2.070e-02 4.444e-02 1.650e+02 0.466 0.642
## SEXmale 2.830e-02 2.724e-02 1.650e+02 1.039 0.300
## Age.SDM8 4.788e-04 1.443e-03 1.650e+02 0.332 0.740
## DX.ABDementia_1 -1.251e-02 3.453e-02 1.650e+02 -0.362 0.717
## DX.ABMCI_1 -3.053e-02 3.113e-02 1.650e+02 -0.981 0.328
## EU.dist 2.993e-03 6.196e-02 1.650e+02 0.048 0.962
##
## Correlation of Fixed Effects:
## (Intr) FC.qnt SEXmal A.SDM8 DX.ABD DX.ABM
## FC.quantile -0.098
## SEXmale -0.009 -0.010
## Age.SDM8 -0.024 -0.034 0.114
## DX.ABDmnt_1 0.000 0.002 0.326 -0.037
## DX.ABMCI_1 0.032 -0.005 -0.035 -0.100 0.549
## EU.dist -0.969 0.104 -0.055 -0.214 -0.040 -0.048
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## `geom_smooth()` using formula = 'y ~ x'
Analysis 5: epicenter connectivity-based prediciton of longitudinal
change rates of SV2A wsocre
# define function
get.epi.pred.beta.followup.shuffled <- function(epicenter.rois, mean.sv2a.w, mean.sv2a.w.cr){
mapply(function(j){
fc.mat.shuffled <- shuffled.fc[[j]]
epi.fc.mat.shuffled <- as.matrix(fc.mat.shuffled[epicenter.rois,])
epi.fc.mat.shuffled[which(epi.fc.mat.shuffled==Inf)] <- NA
mean.epi.fc.mat.shuffled <- colMeans(epi.fc.mat.shuffled,na.rm = T)
df.test.ad.shuffled <- data.frame(epi.fc = mean.epi.fc.mat.shuffled,
sv2a.wscore=mean.sv2a.w,
sv2a.wscore.cr = mean.sv2a.w.cr)
mod.lm.shuffled <- lm(data = df.test.ad.shuffled,sv2a.wscore~epi.fc+mean.sv2a.w.cr)
mod.lm.shuffled <- summary(lm.beta(mod.lm.shuffled))
beta.val <- mod.lm.shuffled$coefficients[2,2]
beta.val
}, 1:1000)
}
# group-level analyses
df.ad.followup <- subset(df.ad.epicenter, !is.na(CR.W.SDM8.V1))
mean.w.bl <- colMeans(df.ad.followup[,paste0("W.SDM8.V",1:200)])
mean.w.cr <- colMeans(df.ad.followup[,paste0("CR.W.SDM8.V",1:200)])
epicenter.bl <- get.epicenter.by.rank(mean.w.bl)
epicenter.bl <- which(epicenter.bl==1)
epi.fc.mat.bl <- fc.mat[epicenter.bl,]
mean.epi.fc.bl <- colMeans(epi.fc.mat.bl,na.rm = T)
beta.shuffled.followup <- get.epi.pred.beta.followup.shuffled(epicenter.bl,mean.w.bl,mean.w.cr)
df.test.followup <- data.frame(epicenter.fc = mean.epi.fc.bl,
sv2a.wscore.cr = mean.w.cr,
sv2a.wscore.bl = mean.w.bl)
mod.lm <- lm(data = df.test.followup,sv2a.wscore.cr~epicenter.fc+sv2a.wscore.bl)
mod.lm <- summary(lm.beta(mod.lm));print(mod.lm)
##
## Call:
## lm(formula = sv2a.wscore.cr ~ epicenter.fc + sv2a.wscore.bl,
## data = df.test.followup)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.279083 -0.066378 0.003909 0.073033 0.245720
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -0.25229 NA 0.00753 -33.505 < 2e-16 ***
## epicenter.fc -0.07310 -0.20377 0.02503 -2.920 0.00391 **
## sv2a.wscore.bl 0.01209 0.02965 0.02845 0.425 0.67139
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1056 on 197 degrees of freedom
## Multiple R-squared: 0.04299, Adjusted R-squared: 0.03327
## F-statistic: 4.425 on 2 and 197 DF, p-value: 0.01319
beta.val <- mod.lm$coefficients["epicenter.fc","Standardized"]
p.val.shuffled <- length(which(beta.shuffled.followup< beta.val))/length(beta.shuffled.followup)
# subject-level analyses
beta.sub.followup <- ddply(df.ad.followup,.(SubID), function(x){
mean.w.sub.cr <- as.vector(as.matrix(x[,paste0("CR.W.SDM8.V",1:200)]))
mean.w.sub.bl <- as.vector(as.matrix(x[,paste0("W.SDM8.V",1:200)]))
epicenters.bl.sub <- which(x[,paste0("Epicenter.bin.V",1:200)]==1)
epi.fc.mat.sub <- fc.mat[epicenters.bl.sub,]
mean.epi.fc.sub <- colMeans(epi.fc.mat.sub,na.rm = T)
df.sub.tmp <- data.frame(mean.w.sub.cr = mean.w.sub.cr,
mean.w.sub.bl = mean.w.sub.bl,
mean.fc.sub = mean.epi.fc.sub)
beta.sub <- coefficients(lm.beta(lm(data=df.sub.tmp,mean.w.sub.cr~mean.fc.sub+mean.w.sub.bl)))[2]
c(beta.sub)
})
beta.sub.followup$variable <- "beta"
t.mod <- t.test(beta.sub.followup$mean.fc.sub)
t.mod
##
## One Sample t-test
##
## data: beta.sub.followup$mean.fc.sub
## t = -1.0008, df = 16, p-value = 0.3318
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.04683032 0.01679366
## sample estimates:
## mean of x
## -0.01501833
Analysis 6: effect of plasma ptau concentration on
connectivity-mediated synaptic loss propagation
# determine subject-specific epicenters for all participants
df.epicenters.subject.all <- ddply(df,.(SubID), function(x){
mean.w.sub <- as.vector(as.matrix(x[,paste0("W.SDM8.V",1:200)]))
epicenters.sub <- get.epicenter.by.rank(mean.w.sub)
epicenters.sub
})
colnames(df.epicenters.subject.all)[2:201] <- paste0("Epicenter.bin.V",1:200)
df.epicenter <- merge(df,df.epicenters.subject.all, by ="SubID")
# derive subject-level connectivity-based synaptic loss
beta.sub.all <- ddply(df.epicenter,.(SubID), function(x){
mean.w.sub <- as.vector(as.matrix(x[,paste0("W.SDM8.V",1:200)]))
epicenters.sub <- which(x[,paste0("Epicenter.bin.V",1:200)]==1)
epi.fc.mat.sub <- fc.mat[epicenters.sub,]
mean.epi.fc.sub <- colMeans(epi.fc.mat.sub,na.rm = T)
df.sub.tmp <- data.frame(mean.w.sub = mean.w.sub,
mean.fc.sub = mean.epi.fc.sub)
beta.sub <- coefficients(lm.beta(lm(data=df.sub.tmp,mean.w.sub~mean.fc.sub)))[2]
c(beta.sub)
})
colnames(beta.sub.all)[2] <- "connectivity.SV2A.beta"
t.test(beta.sub.all$connectivity.SV2A.beta)
##
## One Sample t-test
##
## data: beta.sub.all$connectivity.SV2A.beta
## t = -0.1161, df = 149, p-value = 0.9077
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.01150906 0.01023165
## sample estimates:
## mean of x
## -0.0006387068
df.all.ptau181 <- merge(df.epicenter, beta.sub.all, by ="SubID")
df.ad.ptau181 <- subset(df.all.ptau181, AV45.bin==1)
# linear regression
mod.all <- lm(data = df.all.ptau181,connectivity.SV2A.beta~ log(pTau181)+ Age.SDM8 +AV45.global.mean+ SEX + DX); summary(lm.beta(mod.all))
##
## Call:
## lm(formula = connectivity.SV2A.beta ~ log(pTau181) + Age.SDM8 +
## AV45.global.mean + SEX + DX, data = df.all.ptau181)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.155639 -0.046858 -0.001261 0.043153 0.168463
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) 0.0247728 NA 0.0493242 0.502 0.616
## log(pTau181) 0.0048930 0.0719712 0.0057353 0.853 0.395
## Age.SDM8 -0.0002539 -0.0338735 0.0006301 -0.403 0.688
## AV45.global.mean 0.0020351 0.0215027 0.0080230 0.254 0.800
## SEXmale 0.0049042 0.0365165 0.0116550 0.421 0.675
## DXDementia -0.0115818 -0.0750130 0.0135697 -0.854 0.395
## DXMCI -0.0145663 -0.0888641 0.0144731 -1.006 0.316
##
## Residual standard error: 0.06822 on 143 degrees of freedom
## Multiple R-squared: 0.01597, Adjusted R-squared: -0.02532
## F-statistic: 0.3867 on 6 and 143 DF, p-value: 0.8866
mod.ad <- lm(data = df.ad.ptau181,connectivity.SV2A.beta~ log(pTau181)+ Age.SDM8 +AV45.global.mean+ SEX + DX); summary(lm.beta(mod.ad))
##
## Call:
## lm(formula = connectivity.SV2A.beta ~ log(pTau181) + Age.SDM8 +
## AV45.global.mean + SEX + DX, data = df.ad.ptau181)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.137337 -0.045001 -0.002558 0.045269 0.165296
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -0.0111801 NA 0.0570300 -0.196 0.845
## log(pTau181) 0.0079342 0.1150652 0.0066298 1.197 0.234
## Age.SDM8 0.0002119 0.0281121 0.0007275 0.291 0.771
## AV45.global.mean 0.0084790 0.0891837 0.0092927 0.912 0.364
## SEXmale 0.0084778 0.0634181 0.0134060 0.632 0.528
## DXDementia -0.0175174 -0.1240537 0.0151090 -1.159 0.249
## DXMCI -0.0184616 -0.1246824 0.0161495 -1.143 0.256
##
## Residual standard error: 0.06756 on 106 degrees of freedom
## Multiple R-squared: 0.03807, Adjusted R-squared: -0.01637
## F-statistic: 0.6993 on 6 and 106 DF, p-value: 0.6508
# interaction effect between ptau181 and epicenter functional connectivity
## test for pooled sample
df.stat.plasma.ptau181.interaction <- subset(df.stat.sub, !is.na(pTau181))
df.stat.plasma.ptau181.interaction$pTau181.bin <- ifelse(df.stat.plasma.ptau181.interaction$pTau181>= median(df.stat.plasma.ptau181.interaction$pTau181, na.rm = T),1,0)
mod.lmer <- lmer(data = df.stat.plasma.ptau181.interaction,
mean.sv2a ~ FC.quantile*log(pTau181) + SEX + Age.SDM8 + DX + AV45.global.mean + EU.dist +(1|SubID))
## boundary (singular) fit: see help('isSingular')
summary(mod.lmer)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean.sv2a ~ FC.quantile * log(pTau181) + SEX + Age.SDM8 + DX +
## AV45.global.mean + EU.dist + (1 | SubID)
## Data: df.stat.plasma.ptau181.interaction
##
## REML criterion at convergence: -309.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0100 -0.6880 0.0005 0.5983 3.0817
##
## Random effects:
## Groups Name Variance Std.Dev.
## SubID (Intercept) 0.00000 0.0000
## Residual 0.02587 0.1608
## Number of obs: 452, groups: SubID, 113
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -3.404e-01 2.831e-01 4.420e+02 -1.203 0.230
## FC.quantile 4.635e-02 4.124e-02 4.420e+02 1.124 0.262
## log(pTau181) 4.062e-03 7.900e-03 4.420e+02 0.514 0.607
## SEXmale 2.954e-02 1.597e-02 4.420e+02 1.850 0.065 .
## Age.SDM8 -5.000e-04 8.713e-04 4.420e+02 -0.574 0.566
## DXDementia -1.037e-03 1.804e-02 4.420e+02 -0.057 0.954
## DXMCI -1.898e-02 1.928e-02 4.420e+02 -0.985 0.325
## AV45.global.mean -1.011e-03 1.106e-02 4.420e+02 -0.091 0.927
## EU.dist 5.008e-02 3.773e-02 4.420e+02 1.327 0.185
## FC.quantile:log(pTau181) 2.473e-02 2.982e-02 4.420e+02 0.829 0.407
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) FC.qnt l(T181 SEXmal A.SDM8 DXDmnt DXMCI AV45.. EU.dst
## FC.quantile -0.039
## log(pTa181) 0.071 0.012
## SEXmale 0.002 0.000 0.131
## Age.SDM8 -0.120 -0.011 -0.017 0.060
## DXDementia 0.063 -0.002 0.000 0.035 -0.065
## DXMCI 0.035 -0.008 -0.028 -0.175 0.049 0.424
## AV45.glbl.m -0.036 0.002 0.044 0.152 -0.057 -0.035 0.075
## EU.dist -0.971 0.042 -0.043 -0.042 -0.110 -0.077 -0.074 0.000
## FC.q:(T181) -0.002 0.741 0.003 0.007 0.004 -0.005 -0.008 0.002 0.002
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
y.range <- range(df.stat.plasma.ptau181.interaction$mean.sv2a)
ggplot(data=df.stat.plasma.ptau181.interaction)+
aes(x=Quantile ,y=mean.sv2a, color=factor(pTau181.bin))+
geom_point(shape=19,alpha=0.5)+
geom_line(alpha=0.15,aes(group=SubID))+
geom_smooth(method = "lm",aes(x=Quantile.num , y = mean.sv2a, color =factor(pTau181.bin)), linewidth=0.5)+
labs(x = "",y = "SV2A-PET W-score") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_rect(color= " black", fill = NA),
axis.text.x = element_text(size=12))+
scale_x_discrete()+
scale_color_nejm(name="Plasma p-tau181",labels=c("<median",">median"))+
guides(fill = "none")
## `geom_smooth()` using formula = 'y ~ x'
## test for abeta positive sample
df.stat.plasma.ptau181.interaction.ad <- subset(df.stat.sub, !is.na(pTau181)& AV45.bin==1)
df.stat.plasma.ptau181.interaction.ad$pTau181.bin <- ifelse(df.stat.plasma.ptau181.interaction.ad$pTau181>= median(df.stat.plasma.ptau181.interaction.ad$pTau181, na.rm = T),1,0)
mod.lmer <- lmer(data = df.stat.plasma.ptau181.interaction.ad,
mean.sv2a ~ FC.quantile*log(pTau181) + SEX + Age.SDM8 + DX +AV45.global.mean + EU.dist +(1|SubID))
## boundary (singular) fit: see help('isSingular')
summary(mod.lmer)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean.sv2a ~ FC.quantile * log(pTau181) + SEX + Age.SDM8 + DX +
## AV45.global.mean + EU.dist + (1 | SubID)
## Data: df.stat.plasma.ptau181.interaction.ad
##
## REML criterion at convergence: -309.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0100 -0.6880 0.0005 0.5983 3.0817
##
## Random effects:
## Groups Name Variance Std.Dev.
## SubID (Intercept) 0.00000 0.0000
## Residual 0.02587 0.1608
## Number of obs: 452, groups: SubID, 113
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -3.404e-01 2.831e-01 4.420e+02 -1.203 0.230
## FC.quantile 4.635e-02 4.124e-02 4.420e+02 1.124 0.262
## log(pTau181) 4.062e-03 7.900e-03 4.420e+02 0.514 0.607
## SEXmale 2.954e-02 1.597e-02 4.420e+02 1.850 0.065 .
## Age.SDM8 -5.000e-04 8.713e-04 4.420e+02 -0.574 0.566
## DXDementia -1.037e-03 1.804e-02 4.420e+02 -0.057 0.954
## DXMCI -1.898e-02 1.928e-02 4.420e+02 -0.985 0.325
## AV45.global.mean -1.011e-03 1.106e-02 4.420e+02 -0.091 0.927
## EU.dist 5.008e-02 3.773e-02 4.420e+02 1.327 0.185
## FC.quantile:log(pTau181) 2.473e-02 2.982e-02 4.420e+02 0.829 0.407
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) FC.qnt l(T181 SEXmal A.SDM8 DXDmnt DXMCI AV45.. EU.dst
## FC.quantile -0.039
## log(pTa181) 0.071 0.012
## SEXmale 0.002 0.000 0.131
## Age.SDM8 -0.120 -0.011 -0.017 0.060
## DXDementia 0.063 -0.002 0.000 0.035 -0.065
## DXMCI 0.035 -0.008 -0.028 -0.175 0.049 0.424
## AV45.glbl.m -0.036 0.002 0.044 0.152 -0.057 -0.035 0.075
## EU.dist -0.971 0.042 -0.043 -0.042 -0.110 -0.077 -0.074 0.000
## FC.q:(T181) -0.002 0.741 0.003 0.007 0.004 -0.005 -0.008 0.002 0.002
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
y.range <- range(df.stat.plasma.ptau181.interaction.ad$mean.sv2a)
ggplot(data=df.stat.plasma.ptau181.interaction.ad)+
aes(x=Quantile ,y=mean.sv2a, color=factor(pTau181.bin))+
geom_point(shape=19,alpha=0.5)+
geom_line(alpha=0.15,aes(group=SubID))+
geom_smooth(method = "lm",aes(x=Quantile.num , y = mean.sv2a, color =factor(pTau181.bin)), linewidth=0.5)+
labs(x = "",y = "SV2A-PET W-score") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_rect(color= " black", fill = NA),
axis.text.x = element_text(size=12))+
scale_x_discrete()+
scale_color_nejm(name="Plasma p-tau181",labels=c("<median",">median"))+
guides(fill = "none")
## `geom_smooth()` using formula = 'y ~ x'